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This paper presents tests of a new method for the simultaneous estimation 
of spacecraft attitude and sensor biases, based on a quaternion estimation 
algorithm minimizing Wahba's loss function. The new method is compared 
with a conventional batch least-squares differential correction algorithm. 

The estimates are based on data from strapdown gyros and star trackers, 
simulated with varying levels of Gaussian noise for both inertially-fixed and 
Earth-pointing reference attitudes. Both algorithms solve for the spacecraft 
attitude and the gyro drift rate biases. They converge to the same estimates 
at the same rate for inertially-fixed attitude, but the new algorithm converges 
more slowly than the differential correction for Earth-pointing attitude. The 
slower convergence of the new method for non-zero attitude rates is 
believed to be due to the use of an inadequate approximation for a partial 
derivative matrix. The new method requires about twice the computational 
effort of the differential correction. Improving the approximation for the 
partial derivative matrix in the new method is expected to improve its 
convergence at the cost of increased computational effort. 
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Introduction 

When determining the three-axis attitude of a spacecraft, it is often necessary to simultaneously 
estimate sensor biases and misalignments. An extended Kalman filter or a batch least-squares 
differential correction procedure is generally used for this process [1], These methods, collectively 
referred to as state estimation methods, require that the nonlinear estimation problem be linearized 
about a priori estimates of the attitude, biases, and misalignments. The purpose of this paper is to 
compare the standard batch least-squares differential correction procedure with a new algorithm 
[2] based on the <?-method for minimizing Wahba’s least-squares loss function [3], The new 
algorithm computes the parameters iteratively, but does not linearize about an a priori attitude 
estimate, so it is expected to be more robust than the usual state estimation methods if the problem 
is highly nonlinear or if initial estimates are p<x>r. 

The development of the new algorithm is presented in detail in reference 2, along with some 
historical background, so it will not be repeated here. The following iterative procedure estimates 
the attitude at time t and the vector x comprising the rn non-attitude parameters. An a priori estimate 
x° of x is used to compute the matrix 

Bit, x) =1 a i 0fr, tp x)b i (x)r/ / (x), (1) 

where the unit vectors r ( - are representations in an inertial reference frame of the directions to some 
observed objects, the b,- are the unit vector representations of the corresponding observations in the 
spacecraft body frame, the a ^ are positive weights, and n is the number of observations. The 3x3 
attitude propagation matrix 0(r, r 0 ; x) is the solution of the differential equation 


d0(r, r 0 ; x)/dr = [co(r, x)x] 0(r, r 0 ; x) 

with initial value 

0(% t(), x) - / = the 3x3 identity matrix, 


( 2 ) 

( 3 ) 


where the column vector c o(t, x) contains the components in the body frame of the spacecraft 
angular velocity relative to inertial space. The matrix [vx| is defined for an arbitrary three-vector v 

b y 


(vx] = 


0 v 3 v 2 

v- 3 0 - v, 

v 2 v i 0 


( 4 ) 


The parameters in x may enter the matrix B(t, x) through the kinematics expressed by 0(r, tp x), 
the observation modeling in b^x), or the reference vector models in r^x). The m matrices dB/dxj, 
and the m(m+l)/2 independent matrices d 2 B/dXjdx k expressing the derivatives of B(t, x) with 
respect to the parameters must also be computed. 

Standard methods [4] are next used to compute the largest eigenvalue ^ mar (x) and 
corresponding normalized eigenvector q opt (t, x) of the symmetric 4x4 matrix 


x) = 


B(t, x) + B T (t, x) - / tr B(t, x) 
P 7 ’(f, x) 


p(f, x) 
tr B(t, x) 


( 5 ) 
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with the three-component column vector p(t, x) defined by 

[p(r, x)x] = B T (t, x) - B(t, x). (6) 

Then the optimal attitude matrix for parameter vector x, A op( (t, x), is computed from q opt (t, x) by 
A opt {t, x) = (q 2 - Q r Q) / + 2QQT- 2q[Qx], (7) 

where the three- vector part Q and scalar part q of the quaternion q^f, x) are given by 

qcp/ r ( f » x ) = ( g ) 

The parameter vector update is given by 

Sx opt = ^- 1 (x)[h(x) - W*>(x - x°)], (9) 

where, for j, k = 1, ... , m, 

Wj k (\ ) = [VT° - N T (t, x)M~\t, x)N(t , x)]j k - tr [A opt T (t, x)d 2 B(t, x)/dxjdx k ], (10) 
hj{x ) = tr{[35(r, x)/dxj]A op( T (t, x)}, (1 1) 

Ny(t, x) = {[ dB(t , x)/dXj]A op T(t, x )} 23 - {[dB(t, x)/dxj]A opt T (t, x)} 32 , (12a) 

A^ 2 y(f, x) = { [ dB(t, x)/dx j}A opl T (t, x)) 31 — {[3fi(t, x^/dxj\A 0 pT(t, x)} j 3 , (12b) 

x) = { [95(f, x)/dxj\A op T(t, x) ) | 2 — [\dB(t,x}/dxj]A 0 pT(t, x)} 2 j, (12c) 

and 

M(t, x) = A mflJC (x) / - B(t, x)A opt T {t, x). (13) 


In these equations x° is the a priori estimate of x and W° is a symmetric positive- semidefinite 
matrix of weights assigned to this estimate; it is permissible to assign zero weights to the a priori 
estimate by setting W° = 0. The update dx opt is added to x to get the new parameter estimate. If the 
update is small enough, the procedure is complete; otherwise the computations are repeated from 
equation (1) until convergence is achieved. 

The attitude covariance Pqq, the parameter covariance P n , and the cross-covariance Pfa of 
the converged estimate can be computed as follows: 


/>ee«) = <r, ol 2 l M-‘«, x) + M-Ml, x)N(t, x)W-'U)NT(t, *)]. (14 

p xx = °to? 0 5 : 

and 

/ , 0 X (O = P x e T (t) = o to 2 M-\t, x)N(t, x)W-l(x), (16; 

where « 

or 2 )' 1 07: 

with <J; 2 equal to the variance of the i^ vector measurement. The covariance computation assumes 


the weights to be 


a i ~ (J to?'l G i 2 
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and 


WO = o to 2(P°)-K (19) 

where is the covariance of the a priori parameter vector estimate. An expression for A/ -1 (r, x) is 
derived in the Appendix, eliminating the need for a numerical matrix inversion. 

Application to Gyro Drift Determination 

For the application to be treated in this paper, we assume that the kinematic information for 
attitude propagation is obtained from three gyros aligned with the spacecraft body axes. In this case 
the estimation algorithm assumes that the body rate vector co(r, x) is 

®(r, x) = Q)g(t) - x, (20) 

where a>g(r) is the column vector of gyro outputs and x, a three-component vector of gyro drifts, 
is the vector of parameters to be estimated. These parameters enter B(t, x) through the attitude 
propagation matrices <I)(r, tf, x). The First and second partial derivatives of <b(r, tf, x) with respect 


to the components of x are needed to evaluate the corresponding partial derivatives of B(t, x). The 
partial derivative of equation (2) with respect to Xj is, using equation (20), 

d[d<J>(r, r 0 ; x)/dxj ]/dr = - [co(r, x)x][90(r, t 0 ; \)/dxj ] + [ey x] <X> (t, t 0 ; x), (21) 

where ej is the unit vector along the spacecraft axis. The solution of this differential equation 

consistent with equations (2) and (3) is 

t 

9<D(r, r 0 ; \)/dx: = / <D (r, r'; x)[e • x] <X>(r', r 0 ; x) d t' . (22) 

t Q 

Using the group property and orthogonality of the attitude propagation matrix, 

0(r', r 0 ; x) = <D(r', t; x)0(r, r 0 ; x) = t'; x)<D(r, r 0 ; x), (23) 

and the relation 

Ofey x] O r = [(<De^ )xj , (24) 

which holds for any proper orthogonal matrix O, gives 

d<P(t, t 0 ; x)/dxj - - r\py(r, r 0 ; x)x]0>(r, r 0 ; x), (25) 

where \j /.(/, tpi x) is the column of the matrix 

t 

^(r, t 0 ; x) s - J d>(t, t'; x) dr' = [xjr^r, r 0 ; x) \y 2 {t, t 0 ; x) \j/ 3 (r, r 0 ; x)]. (26) 


The matrices 0(r, rp; x) and ^(r, tp; x) are also used for the usual state estimation methods 
[5, 6]. They are computed by adding up contributions over time intervals r t+1 - 1 ,• , which are 
chosen to be short enough that variations in o) over the interval can be neglected. Thus 

0(r /+1 , tQ, x) = <B(r f+1 , ; x)<D(r / , r 0 ; x), (27) 

where 

0(r /+1 , r,- ; x) = /- [ux]sin|<t>| + Tux] 2 (l — cos|<>|), (28) 

u is the unit vector 

u = co/|co|, (29) 
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and |<|>| is the length of the vector 


(30) 


Similarly, 


4> = ®(r i+ 1 - fy ). 


x) = 'F(t J+1 , fy ; x) + <D(fy +1 , fy ; x^ry, f 0 ; x), 

with 

^(fy+l, *y ; x) = - (t i+l - fy ){/- [uxJI^-^l -cos|<J>|) + [ux]2(l - |<J>|- 1 sin|(t)|)}. 

The second partial derivatives of 0(r, to; x) only appear in VF(x) and not in h(x), so 
approximate forms will be used for these partials. They are also computed by adding up 
contributions over short time intervals fy +1 — fy , giving 


(31) 

(32) 


0 2 O(r /+1> r 0 ; x)/dxj dx k = [0 2 O(fy +1 , fy ; x)/d Xj dx k ]<b(fy, r 0 ; x) 

+ [90(fy +1 , t t ; x)/9xy][3<D(fy, t 0 ; x)/dx k ] + [30(f J+1 , ^ ; x)/dx k ]\d<&(t i , t 0 ; x)/dxj\ 

+ <&(fy +1 , ^ ; x)[a 2 <D(fj, t 0 ; x)/dxj dx k ], (33) 

where the first partial derivatives are give by equations (25) - (32) and where 

3 2 d>(fy +1 , t t ; x)/dxj dx k = (fy +1 - ti ) 2 {\ {tj e k T + e k e?- IbjjJ) 

+ J (4>jle k x] + <p k [ej x] + 8 jk [<\>x]) } . (34) 

The approximation is in equation (34), which is valid to first order in <|>. Starting the iterative 
computations of equations (27), (31), and (33) requires initial values for the matrices: the identity 
for O(f 0 , tQ; x) from equation (3), and zero for ^(tg, f 0 ; x) and d 2 0(rg, r 0 ; x)/dxj dx k . 

Observation Modeling 

Star tracker data are used to estimate the spacecraft attitude and gyro drifts. Each star tracker 
measurement is a two-component vector y y giving the location of the star image in the focal plane 
of the sensor. For attitude estimation with the new method we need to compute the star unit vector 
in the spacecraft body frame by in terms of the measurement data yy. The star unit vector in the 
sensor frame is 

8/ = (1 + lyy l 2 )- 1/2 [(yy) P (y/) 2 , 1] J , (35) 

and then by is given by 

by = CyTsy , (36) 


where Cy is the proper orthogonal 3x3 matrix defining the orientation in the body frame of the star 
tracker making this observation. 

Data Simulation 

Simulated gyro data and star tracker data are used to test the algorithm. The simulation assumes 
a constant angular velocity vector (O true . The gyro data are simulated by adding varying levels of 
Gaussian noise to the components of <% ue . A true attitude matrix is computed by integrating the 
angular rates; 

tiA^yg(t)/d/ = — [WfrufX] A^yg(r), (37) 


with some specified initial attitude matrix A lrue (tQ). 
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A star is initially simulated for each star tracker by randomly generating a measurement vector 
y,- within the star tracker field of view. Equations (35) and (36) then give the star unit vector bj- in 
the body frame, and the star unit vector in the inertial reference frame is given by 

*7 = > (38) 

where fj is the simulation time. For successive simulation times f,- , the reference vector r f - is held 
fixed and the vector in the body frame is computed as 

= A true (ti)r i . (39) 

Then the corresponding vector in the star tracker reference frame is given by the inverse of 
equation (36), 

s / = C,b,- ; (40) 

the measurement vector y,* by the inverse of equation (35), 

y j = (5/) 3 _1 [(^) 1 . (Sj) 2 ] r » (41) 

and Gaussian noise is added to the two components of y ,• . This process is continued until the star 
has been tracked for some fixed number of observations or until it leaves the field of view, at 
which time a new star is randomly placed in the field of view. Earth and Sun interference are 
neglected in these simulations. 

Comparison Algorithm 

The algorithm chosen for comparison is a batch least-squares differential correction algorithm 
similar to that employed in the attitude ground support system of the Upper Atmosphere Research 
Satellite (UARS) [6]. The algorithm provides a least-squares estimate of a six-component state 
vector 

8XT= [80?' 5x7}, (42 ) 

where 50 is the attitude error vector at epoch, and 8x is the error in the gyro drift estimates. This is 
updated iteratively as follows. At the start of each iteration, an estimate x of the gyro drifts and 
A est (tQ) of the epoch attitude are available. For each measurement y,-, a predicted value g,- is 
computed from the known reference vector i*j- by equations identical to (39) - (41), but with the 
unknown attitude matrix A^itj) replaced by 

^est^O = *0» x)A es i(t g). (43) 


The computed value g,- is seen to depend on both x and A^tg). The optimal state update is the 
solution of 


F8X= S c r 2 G i T(y i -g i )- 


0 

(/>0)-l(x - x°) ’ 


where 


n „ „ To 0 1 

F = E <T; + nil’ 

1=1 1 1 1 [o (FVJ 
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with Gf 2 , P°, and x° as defined previously, and with 0 a three-vector of zeros and 0 a 3x3 matrix 
of zeros. The 2x6 matrix G,- of partial derivatives of the errors of the i™ measurement with respect 
to 5X is given by 



- 1 - Oi), 2 

(yih 

1+ O’* 2 


-Wi m 


■i [<&(*,-, to ; *) ^(f,-, ^ x)]. 
This state update gives new estimates of the gyro drifts and attitude: 


(46) 


X/i gw — x + 8x, (47) 

and 

A new (t 0 ) = {/ - iserHsexisinisei + isertsexpa _ cosisei) M^(r 0 )- (48) 


This iterative procedure is repeated until convergence is achieved. An estimate of the covariance 
matrix is provided by 


The initial attitude to begin the first iteration is provided by the q-method, as embodied in equations 
(1) - (7). 

Numerical Examples 

Tests were performed for both inertially-fixed and earth-pointing spacecraft attitudes, with star 
tracker orientations and other parameters corresponding to the Gamma Ray Observatory (GRO) [7] 
and UARS [6] spacecraft, respectively. Two star trackers were modeled with 8 degree by 8 degree 
fields of view and with an angle of approximately 73 degrees between their boresights. Some tests 
were performed with perfect star tracker measurements, but the results presented in this paper all 
include Gaussian noise on each star tracker output with standard deviation of 8 arc seconds, or 
3.88 x 10 -5 radians. The time interval between star tracker measurements was 32.768 seconds, the 
interval used by the UARS onboard computer. The data were simulated with no gyro bias, and the 
estimations were performed with non-zero initial bias estimates, so the bias estimate is the same as 
the bias error for these tests. The initial bias error was 10 -4 radians/sec along either the spacecraft 
roll or yaw axis, but only representative results with initial roll bias errors are given below. 

The epoch time tQ for the estimation was taken to be the time of the first observation. In all but 
four tests, the true attitude matrix at epoch was set equal to the identity matrix. The tests for one 
simulation case were repeated with four different true attitude matrices at epoch: 


^ee p ex =F _i 
^xG p xx_ 


Atru*(t6) ~ 


0.352 0.864 

-0.864 0.152 

0.360 -0.480 


0.360 
0.480 , 

0.800 


^lrue(to) - diag[l, —1, 1], 


^tr««(^) = dia g[- 1 » 1. —11. 


(50a) 


(50b) 

(50c) 
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and 


A truei t o) = diag[-l, - 1 , 1 ], (50d) 

where diag[. . .] denotes a 3x3 matrix with the given elements on the main diagonal and zeros 
elsewhere. The attitude and bias errors for these different initial attitudes were identical to those for 
A true( t o) = I within the precision of the output, as they should be. The covariance P° of the a priori 
bias estimates was taken to be infinite for all tests. 

A representative subset of the tests is presented in Tables 1-7. At least 10 iterations were 
performed in each case, and the errors for all the iterations after those presented are identical to the 
errors of the last iteration in the table, to the precision of the table. The first "iteration" in the 
differential correction (DC) columns is not really a DC iteration; it is an initial attitude estimation 
using the < 7 -method, as explained above. Thus the bias error after one DC "iteration" is the a priori 
error. The last line in each table is the estimate of the error standard deviations from the covariance 
matrices of equations (14) and (15) or equation (49). 

Tables 1-6 present the tests with inertially-fixed attitude. These are in pairs: Tables 1 and 2 give 
the results for the highest observability case with two star trackers and a full orbit of data, Tables 3 
and 4 have two star trackers but only 10 observations in each, and Tables 5 and 6 are for the case 
of only 10 observations in a single star tracker. Each simulated star was observed five times in 
these tests, so the cases in Tables 5 and 6 contain only two stars; the angular separation between 
these stars was 1.3 degrees. In each pair of tables, the first (odd-numbered) presents the results 
with no gyro noise, and the second (even-numbered) shows the effects of Gaussian noise on each 
gyro with standard deviation of 1 degree/hour, or 4.848 x 10-6 rad/sec. 

The most important aspect of the tests, as concerns this paper, is the comparison of the results 
of the new method to those of the DC. The bias and attitude errors are not the same at each 
iteration, but both the general rate of convergence and the final converged errors are almost 
identical. Where there are differences, the errors of the new method are slightly lower, but not by a 
significant amount. 

In the cases without gyro noise, the covariance matrix is a good indicator of the estimation 
errors. This correspondence is especially striking in Tables 1 and 3, while the actual errors in Table 
5 are about 20 times less than the covariance matrix would indicate. The errors in the latter case are 
remarkably small considering the poor measurement geometry, with only two reference vectors 
separated by 1.3 degrees. When gyro noise is included, the actual errors can exceed the covariance 
estimates; this is not surprising since the covariance computation does not take gyro errors into 
account, nor does any other part of the estimation process. When unrealistically large gyro noise 
with standard deviation of 100 degree/hour was included, both estimation procedures became 
unreliable. The new method failed catastrophically when the nominally positive-definite matrix W 
defined by equation (10) developed a negative element on its main diagonal. The DC did not 
become singular, since the matrix F of equation (45), unlike W, is manifestly positive-semidefinite; 
but the bias estimation error increased monotonically for the first 10 iterations. Thus the new meth- 
od is somewhat less robust than the DC in this case; but this is not very significant since a Kalman 
filter or smoother should probably be used in the presence of large amounts of dynamic noise [5]. 
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Table 1. Bias and Attitude Errors for Inertially-fixed Attitude 
with Two Star Trackers, 95.6 Minutes of Data, and no Gyro Noise 



Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

1.00D-4 

2.88D-1 

1.68D-6 

2.88D-1 

2 

5.65D-7 

1.50D-3 

1.69D-9 

4.83D-3 

3 

1.68D-9 

5.83D-6 

1.68D-9 

5.82D-6 

4 

1.68D-9 

5.81D-6 

1.68D-9 

5.80D-6 

Covariance 

2.90D-9 

9.53D-6 

2.91D-9 

9.54D-6 



Table 2. Bias and Attitude Errors for Inertially-fixed Attitude 


with Two Star Trackers, 95.6 Minutes of Data, and Gyro Noise of 1 deg/hour 


Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

1.00D— 4 

2.87D-1 

1.73D-6 

2.87D-1 

2 

6.54D-7 

2.66D-3 

2.90D-7 

5.46D-3 

3 

2.90D-7 

1.98D-3 

2.90D-7 

1.98D-3 

Covariance 

2.90D-9 

9.53D-6 

2.91D-9 

9.54D-6 


Table 3. Bias and Attitude Errors for Inertially-fixed Attitude 



with Two Star Trackers, 5.5 Minutes of Data, and no Gyro Noise 


Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

l.OOD^l 

1.48D-2 

2.12D-7 

1.48D-2 

2 

2.17D-7 

2.95D-5 

2.07D-7 

3.02D-5 

3 

2.07D-7 

2.93D-5 

2.07D-7 

2.93D-5 

Covariance 

2.10D-7 

3.67D-5 

2.10D-7 

3.68D-5 


Table 4. Bias and Attitude Errors for Inertially-fixed Attitude 
with Two Star Trackers, 5.5 Minutes of Data, and Gyro Noise of 1 deg/hour 


Batch DC 

New Method 

Iteration Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 1.00D— 4 

1.49D-2 

3.50D-6 

1.49D-2 

2 3.50D-6 

2.71D-4 

3.50D-6 

2.71D-4 

3 3.50D-6 

2.72D-4 

3.50D-6 

2.72D^1 

Covariance 2.10D-7 

3.67D-5 

2.10D-7 

3.68D-5 
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Table 5. Bias and Attitude Errors for Inertially-fixed Attitude 
with One Star Tracker, 5.5 Minutes of Data, and no Gyro Noise 



Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

1.00D-4 

1.63D-2 

1.14D—4 

1.63D-2 

2 

2.29D-5 

3.25D-3 

1.07D-6 

5.69D-3 

3 

1.38D-6 

2.94D-4 

1.08D-6 

2.49D-4 

4 

1.10D-6 

2.55D— 4 

1.08D-6 

2.51D-4 

5 

1.11D-6 

2.55D-4 

1.08D-6 

2.51D-4 

Covariance 

2.39D-5 

4.18D-3 

2.40D-5 

4.19D-3 


Table 6. Bias and Attitude Errors for Inertially-fixed Attitude 



with One Star Tracker, 5.5 Minutes of Data, and Gyro Noise of 1 deg/hour 


Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

1.00D-4 

1.44D-2 

1.18D-4 

1.44D-2 

2 

2.25D-5 

1.16D-2 

4.39D-6 

7.11D-4 

3 

4.51D-6 

1.54D-2 

4.41D-6 

1.54D-2 

4 

4.50D-6 

1.54D-2 

4.41D-6 

1.54D-2 

Covariance 

2.37D-5 

4.16D-3 

2.38D-5 

4.18D-3 


Table 7. Bias and Attitude Errors for Earth-Pointing Attitude 



with One Star Tracker, 5.5 Minutes of Data, and no Gyro Noise 


Batch DC 

New Method 

Iteration 

Bias (rad/sec) 

Attitude (rad) 

Bias (rad/sec) 

Attitude (rad) 

1 

1.00D— 4 

3.01D-2 

8.04D-5 

3.01D-2 

2 

1.37D-5 

2.54D-3 

1.73D-5 

3.02D-2 

3 

4.74D-6 

8.52D-4 

7.17D-6 

3.05D-3 

4 

4.88D-6 

8.76D-4 

5.31D-6 

1.27D-3 

5 

4.88D-6 

8.76D-4 

4.97D-6 

9.48D^1 

6 

4.88D-6 

8.76D-4 

4.91D-6 

8.90D— 4 

7 

4.88D-6 

8.76D-4 

4.89D-6 

8.79D-4 

8 

4.88D-6 

8.76D-4 

4.89D-6 

8.78D^1 

9 

4.88D-6 

8.76D-4 

4.89D-6 

8.77D-4 

Covariance 

1.24D-5 

2.15D-3 

1.12D-5 

1.96D-3 
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The attitude errors in Table 2 are larger than those in Table 4, showing that estimators that do not 
handle dynamic noise correcctly should not be used with long spans of data including the effects of 
dynamic noise. 

The tests with Earth-pointing attitude used a constant pitch rate of -1.083073 x lO - ^ rad/sec, 
which corresponds to an orbit period of 96.7 minutes. For these tests, a simulated star was tracked 
until it left the star tracker field of view. The test parameters were otherwise the same as for 
inertially-fixed attitude. The new method did not fare as well in these tests; it generally required 
more iterations than the DC to converge, although the final errors were virtually identical. This 
suggests the presence of errors in the matrix W(x) that steers the estimates to their optimal values, 
and not in the vector h(x) that identifies the optimum once it has been reached. Table 7 presents the 
results of a test with a single star tracker, observing only two stars with angular separation of 7 
degrees. This is a particularly discouraging example, in which the DC converged in four iterations, 
while the new method required nine. Since this is a low observability case, in which attitude 
kinematics information is more important compared to the measurements than in a high 
observability case, an accurate computation of VT(x) is especially important. 

The greater success of the new method for inertially-fixed attitude than for non-inertial attitude 
suggests the inadequacy of the approximation of equation (34) for the matrix of second partial 
derivatives 3 2 <E>(r 1 * +1 , tf, x)/dx '• dx^, which appears in lf(x) and not in h(x). This approximation 
should be replaced by one that is valid for all values of the rotation angle <J), subject to the 
assumption that the angular rates are approximately constant between observations. This may also 
avoid the failure of the new method in the test with 100 degree/hour gyro noise, since this has an 
effect on the propagation of the partial derivative matrices similar to the effects of actual angular 
rates of the same size. 

The computational effort required by the two algorithms was also measured. Both algorithms 
were implemented in double-precision Fortran and executed on a DEC VAX 1 1/780. The CPU 
times were proportional to the number of iterations performed, the times per iteration for the two 
methods being 

t nDT , =13 + 15.2 n msec (51a) 

CPU , new ' 7 

and 

l CPU DC = +6.7 n msec. (51b) 


where n is the number of observations. The coefficient of n in these times can be interpreted as the 
time required to process a measurement, including propagation of the attitude transition matrix, 
partial dervative matrices, and so forth. The ^-independent term represents the end-of-iteration 
computations, including matrix inversions and computation of updates to the bias vector and 
attitude matrix. Thus equation (51) shows that the measurement processing is more expensive for 
the new method, while the end-of-iteration computations of the DC require more effort. The exact 
CPU times will vary from case to case, but the DC appears to be about twice as fast, overall, as the 
new method for the numbers of measurements typically processed. Improving the computation of 
the matrix of second partial derivatives for the new method will require even more effort to process 
each measurement. 
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Conclusions 

These tests establish the validity of a new method for the simultaneous estimation of spacecraft 
attitude and sensor biases, based on a quaternion estimation algorithm minimizin g Wahba's loss 
function. The new algorithm performs as well as a batch least-squares differential correction in 
tests with inertially-fixed attitude, in the sense of converging to equally accurate estimates in the 
same number of iterations. The new algorithm converges more slowly than the differential 
correction for Earth-pointing attitude, probably owing to the use of an inadequate approximation 
for a partial derivative matrix in the new method. The new method does not show any advantages 
in terms of robustness or speed of convergence, and in addition requires about twice the 
computational effort of the differential correction. It is hoped that improving the approximation for 
the partial derivative matrix in the new method will improve its convergence and/or robustness, 
without adding significantly to its computational burden. 

Appendix 

The matrix B(t, x) has the singular value decomposition [8] 

B = U + S'VJ, (Al) 

where U + and V + are proper orthogonal matrices, and 

S' = diag[5 1 ,5 2 , 5 3 ], (A2) 

a 3x3 matrix with Sj, S 2 , and S 3 on the main diagonal and zeros elsewhere. The arguments t and 
x are omitted from this and all subsequent equations for notational simplicity. The optimal attitude 
estimate is given in terms of these matrices by [8] 

^opt — U+V+T- (A3) 


The maximum eigenvalue of the matrix K defined by equation (5) is related to the optimal 
attitude by [2] 


^max 

(A4) 

where tr denotes the trace. Equations (Al) - (A4) and (13) show that 


detfl =S 1 S 2 S 3 , 

(A5) 

^max ~ S l +S 2 + ^3’ 

(A6) 

M = U + diag[S 2 + S 3 , S 3 + S v S x + S 2 ]Uj. 
and 

(A7) 

det M = (S 2 + S 3 )(S 3 + S l )(S l + S 2 ). 

(A8) 

We now define the scalar 


A little algebra shows that 

(A9) 

k = S 2 S 3 + S 3 Sj +SjS 2 , 

(A10) 
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and 


k X max - det B = det M, 


(All) 


Kl + BBT = u + diag[(S 3 + 5 1 )(5 1 + S 2 ), (5! + S 2 )(S 2 + S 3 ), (S 2 + S 3 )(S 3 + = adj M, 

(A12) 

where adj denotes the adjoint matrix. Equations (A1 1) and (A12) give the desired result 

A/- 1 = (k h max - det B)~\k I + BBT). (A13) 

The evaluation of A/ -1 by means of equations (A9) and (A13) does not require the singular value 

decomposition of B. 
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